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I. INTRODUCTION 


A. BACKGROUND 

A primary goal of a computer graphics system is to provide the user with 
different views of objects. Sometimes the objects that the user manipulates are 
simple in nature and can be constructed easily with the primitives provided by 
the graphics support package. However. in most applications areas the user is 
concerned with more complex objects. The display of three-dimensional surfaces is 
one such application area in computer graphics and is the area that this study 
explores. 

It is the primary intent of this study to stimulate the reader’s interest in the 
area of three-dimensional surface generation and display. To provide this 
stimulation. we combine the power of certain mathematical techniques and a high 
performance graphics environment to design and implement a set of functions 
that can be used to create, manipulate, and display three-dimensional solid-filled 
surfaces. Once developed, the reader will not only be able to use these functions 
to explore the design, representation, and rendering of such surfaces but also will 
be able to use these functions in other fields that can derive benefit from their use 


such as cartography, robotics, computer vision, and artificial intelligence. 


B. PROJECT DESCRIPTION 
1. Motivation 

For developing graphics applications. the typical graphics environment 
provides the user with a variety of sophisticated, powerful tools and a high-level 
language. Almost inevitably, an application calls for something that is not 
provided for or if provided is not acceptable. The lack of solid-filled curved 
surface support in one particular high-performance graphics workstation 
prompted this study. By providing this capability we can create applications that 
enhance and expand what we are able to do and conceive of doing on these 
workstations. 

2. Proposed Capabilities 

To provide this capability for generating solid-filled parametric bicubic 
surface patches, we must look for a method that is simple, powerful, 
understandable, and implementable in software. One approach would be to start 
from scratch and develop a solution. This, however, is usually not a wise way to 
proceed. A solution derived in this manner is likely to have severe limitations 
such as being computationally inefficient, lacking robustness, and difficult to use. 
Another more prudent way to proceed is to find some basic mathematical 
techniques that can be applied to the problem in a form consistent with the 
existing graphics environment. Solutions developed in this manner are more 
easily accepted since their basis lies in proven mathematical techniques and 


generally meets the requirements we seek of being simple, powerful, 
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understandable and implementable. We chose to develop our new capability 
following the latter method - using pre-existing mathematics. 
a. Parametric Bicubic Surface Construction 
One of the primary capabilities needed when working with surfaces is 
to have an effective method for describing, manipulating and displaying them. 
For this, we base our work with surfaces on a mathematical model - the 
parametric bicubic surface patch. Choice of this method allows the user to rapidly 
develop a smooth surface and display it as a wireframe image through the 
specification of only a few points called control points. 
b. Polygonal Parametric Bicubic Surface Decomposition 
Another capability required is to extract, from the mathematical 
model, the information needed to construct a solid surface. For this, we need to be 
able to decompose any surface represented in parametric bicubic form into 
arbitrarily small polygons. These polygons can then be used in the construction 
of a polygon mesh. This capability permits the user to use standard or 


customized algorithms to manipulate the surface. 


C. PROGRAMMING ENVIRONMENT 

The IRIS Turbo 2400 Graphics Workstation, manufactured by Silicon 
Graphics, Inc., is the target programming environment for the design, 
development, and implementation of the parametric bicubic surface constructor 


and polygonal bicubic surface decomposition functions. The IRIS has special- 


purpose hardware that is designed to replace less efficient software. The system 
supports real-time color graphics, the Unix operating system software. and the C 
programming language. A high resolution color monitor provides the output 
device for displaying the graphical output of the modeling and user defined 
functions. 

In addition to the standard programming support provided by the UNIX 
operating system, the IRIS has an extensive collection of utility and graphics 
functions contained in a Graphics Library. This library provides the user high- 
level access to the hardware, enabling graphical objects to be easily manipulated 
as geometrical objects (points, lines, polygons, etc.) rather than pixels. A series of 
coordinate systems and mapping instructions also provides the user with the 


capability to define such objects in world coordinate space. 
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Il. THREE DIMENSIONAL REPRESENTATION OF SURFACES 


We stated in the previous chapter that an important application area in 
computer graphics is concerned with the three-dimensional representation of 
surfaces but we did not describe how one might do this. There are many ways to 
represent surfaces. The two most commonly used representations are: polygon 


meshes and parametric bicubic patches. 


A. POLYGON MESHES 

A polygon mesh is nothing more than a set of connected polygonally planar 
surfaces. There are several ways in which these planar polygons can be 
represented - explicit polygons, vertex list pointers, explicit edges, etc.. Each of 
these representations has its own advantages and disadvantages and various 
criteria can be applied to evaluate the representation. Criteria such as how much 
primary and secondary storage is available, how easy is it to identify the polygons 
sharing an edge, how difficult is it to display the mesh, name only a few 
commonly used for evaluation. Regardless of how the mesh is represented, there 
are many algorithms available for processing it. For example, algorithms have 
been developed to remove hidden surfaces from an object, while others can 


produce lighting and shading effects on the surface of the object. 
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Many objects that have planar features such as buildings, tables, and cabinets 
can easily be represented by a polygon mesh. However, polygon meshes are not 
limited to representing only planar objects. They can also be used to represent 
curved surfaces. To represent a curved surface with a polygon mesh, one creates 
an approximation to the surface by using arbitrarily small polygons. The smaller 
the polygons the better the approximation. However, certain difficulties arise 
when representing surfaces as polygon meshes. As one approximates the surface 
with smaller and smaller polygons, both space and execution time of the 


algorithms that process the mesh increase linearly. 


B. PARAMETRIC BICUBIC PATCHES 

Parametric bicubte patches are mathematical models of a surface and 
represent one of the simplest mathematical elements we can use to model an 
arbitrary surface. A patch can be defined as a curve bounded collection of points 
whose coordinates are given by a continuous, two-parameter, single valued 


function in the form 


x = x(s,t) 
y = y(s,t) 
z = z(s,t) 


where the parameters are restricted to 


s,t « [0,1]. 
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Patches are used to describe the surface of an object in a piecewise manner 
much like the polygon mesh. That is, many patches are used to describe the 
surface of an object and the total surface is generated by displaying all of its 
individual patches. 

Bicubic patches are most often used to generate line drawings of three- 
dimensional objects - often called wireframe representations. One benefit of using 
bicubic patches is that fewer bicubic patches than polygonal patches are needed 
to represent a curved surface to a given accuracy. On the other hand, the 
algorithms for working with bicubic surfaces are more complex than those for 
polygon meshes. In addition, wireframe representations can exhibit certain 
deficiencies. For example, objects that are represented by wireframes are often 
ambiguous. That is, you can look at the same object and get different visual 
interpretations. One reason for the different interpretations is the ability to see 
through a wireframe representation of an object. While not a problem for some 


applications, it can be a serious problem for others. 
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lil. PARAMETRIC CUBIC CURVES AND SURFACES 


A. GENERAL 
The method that we use to generate surfaces is based on using parametric 
cubic curves so it is helpful to review the mathematical basis for these curves. A 
parametric cubic curve has the property that x, y, and z can be defined as third- 
order polynomials for some variable t [Ref. 1:pp. 34] : 
x(t) = a,t® + b,t? + ¢,t + d,, 
V(Ujeauat? tab to 4 cl tide, 
z(t) = a,t® + b,t? +. ¢,t + d,. 
These equations are known as the algebrate form of a parametric cubic curve. 
In this form, we can identify a unique set of 12 constant algebraic coefficients. 
These algebraic coefficients determine a unique parametric cubic curve; they 
determine the size and shape of the curve and its position in three-dimensional 
space. Two curves of the same shape have different algebraic coefficients if they 
occupy different positions in three-dimensional space. Because we want to deal 
with a finite segments of the curve, we limit the range of the parameter, without 
loss of generality, to O0<t<1. We call these finite pieces curve segments. A curve 
segment is nothing more than a point-bounded collection of points. In our case 
the points are bounded at t=O and t=1. Each equation in the algebraic form can 


be expressed as a vector product as follows: 
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Using the vector product form is usually more convenient as it separates the 
distinct parameters of the parametric equation into unknown coefficients of x(t) 
and the parameter t that we wish to manipulate. Here T is the row vector of 
powers of t, while C, is the column vector of coefficients of x(t). Similarly the 
parametric equations for y(t) and z(t) can be written as y(t) = TC, and z(t) = TC,. 
By varying the parameter t from 0 to 1 in each equation we define the curve 
segment. 

Arbitrarily assigning values to these unknown coefficients results in defining a 
curve in three-dimensional space. However, it is not easy to determine the 
properties of this curve. What we wish to do is establish some constraints on 
these coefficients. We want the curves we generate to have some predictable 
properties. To solve the equations for these unknown algebraic coefficients, we 
establish a set of constraints, thereby defining a unique cubic curve with 
predictable properties. To illustrate this process we look at some example cubic 


curves for which the constraints are well-known. 


B. CUBIC CURVE EXAMPLES 
1. Hermite Curve 
The Hermite cubte curve is determined from its endpoints (P,, P,) and 


endpoint tangents (R,, R,). In the literature [Ref. 2:pp. 516-519} (Ref. 3:pp. 123- 
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129], we find the geometric form of the Hermite curve to be 
Qnty = TMiGn- 


In geometric form, T is the row vector of the powers of the parametric 
variable t, M, is the basis matrix, and G, is the geometry vector. A basis matrix 
refers to a constraint procedure that is embodied in matrix form and a geometry 
vector contains the control points used to guide the curve. 


In this particular case. the Hermite basis matriz (M,,) is 


2 = 2h ee 
Se 
Oo 0 Fs 
1 O O O 


and the Hermite geometry vector is 


P, 
P, 
R, 
R, 

Now using the above formulation, given two points and their tangents, we 
can evaluate x(t), y(t), and z(t) for O0<t<1 and find all points on the Hermite form 
of the cubic curve from P, to P, with starting tangent vector R, and ending 
tangent vector R,. It is through these constraints (M,) that the control points 


(G,) control the parametric equations and produce an equation that can generate 


a discretely sampled curve segment in three-dimensional space. 
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If we take the product TM,,. we have 
TM, = | (2t°—3e?44) (—2t*+st?) (t°—2t7+t) (t®—t?) |. 


These four functions of t in the product TM, are often called blending functions 
[Ref l:pp. 48-52]. As the name implies, they blend the effects or contributions of 
the endpoints and tangent vectors to produce the intermediate point coordinate 
values over the domain of t. 
2. Bezier Curve 

The defining form for a Bezter cubic curve is similar to the Hermite form. 
The difference is in the definition of the endpoint tangent vectors. The Bezier 
form uses four points ( P,, Pz, Ps, P, ) instead of 2 points and 2 tangent vectors. 
The tangent vectors at the endpoints in Bezier form are determined by the line 
segments P,P, and P,P, . The Bezier cubic curve passes through the first and 
fourth control points (P, and P,) and uses the second and third points (P, and 
P,) to determine the shape of the curve. [Ref. 1:pp. 113-125] [Ref. 2:pp. 519-521]. 


The geometric form of the Bezier curve is 
Q,(t) = TM,G, 


where the Bezter basis matriz (M,) is 


eo an 

3 -6 3 O 

M, —3 3 O O 
1 O O O 
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and the Bezier geometry vector is 


The Bezier form of the cubic curve is more widely used in computer 
graphics than the Hermite form. A primary reason for its popularity is that the 
geometry matrix of four points (G,) is more intuitive for an interactive user. The 
user has only to manipulate the four points and does not have to specify the 
tangent vectors. It is usually easier for a person to think about manipulating 
points rather than trying to manipulate points and tangent vectors. 

3. Other Useful Cubic Curves 

The Hermite and Bezier cubic curves are not the only forms of cubic 
curve that are available. Two others are the Cardinal Spline and the B-Spline. 

a. Cardinal Spline 

The Cardinal Spline curve passes through the two interior control 
points (P, and P;) and uses the points P, and P, to define the shape of the curve 


[Ref. 4:p. 11-4]. The geometric form of the Cardinal curve is 
Q.(t) m I'M,G, 


where the Cardinal basis matriz (M,) is 
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—a 2-a a-2 a 
Za a-—3 3—2a -—a 
—a 0 a 0 

0 1 0 0 


and the Cardinal geometry vector is 


P, 
P, 
Ps; 
P, 


The scalar coefficient a in the Cardinal basis matrix must be positive and 
determines the length of the tangent vector at point P. and Ps. 
b. B-Spline 


The geometric form of the B-Spline curve is 


Q,,(t) = TM,.Gos 


where the B-Spline basis matriz (M,,) is 


—i) 3 —3 yl 
3 —-6 3 O 
-—-3 O 3 O 
1 4 1 #O 


and the B-Spline geometry vector is 


P, 
Ps 
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In general, the B-Spline curve does not pass through any control 
points but is continuous and also has continuity of tangent vectors and of 
curvature (that is, first and second derivatives are continuous at the endpoints). 
The Hermite and Bezier forms have only first-derivative continuity at the 
endpoints, but do pass through control points [Ref. l:pp. 125-146] [Ref. 2:pp. 


521-523) 


C. DEFINING SURFACES 

By adding a new parameter s and additional algebraic coefficients to the cubic 
curves in the previous section, we can define the algebraic form of a bicubic 
surface patch [Ref. l:pp. 156] as 


i=3 jJ=8 


p(s,t) = © Vays't? 


1=0 J=0 
with the restriction on the parametric variables to 
s,t € [0,1]. 


By varying both parameters from 0 to 1 in each equation, we define al] 
points on the surface patch. Assigning one parameter a constant value and 
varying the other, results in a cubic curve. 

Expanding the above equation in terms of x(s,t) and noting that the terms for 


y(s,t) and z(s,t) are similar we have 
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x(s,t) = as.s*t* + a,.5°t? + ag,5°t + agoS® + agg5*t® + a2.5°t? + a5,87t + aos” 


+ A,,St° + a,ost? + ay St + a,9S + Agst” + anot? + apt + agg- 
Written in vector product form 
Misc) = SC. 1” 


where 


Ass As2 As, Ago 
423 422 421 420 
Arg 412 411 410 


4o3 4292 401 4200 


and T? is the transpose of the matrix T . 


From these equations we can see that there are 48 degrees of freedom or 


algebraic coefficients that we must specify. Like the cubic curve, a change in any 


one of these coefficients defines a different surface. 


The complete algebraic manipulation of the equations to arrive at the 


following equation is similar to that of the curve process described in the previous 


section. For the Bezier surface patch, the geometric form of the equation is: 


x(s,t) = SM,Q,MITT 


where M, is the same M, as in the Bezier curve equation, M{# is its transpose, and 
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Q, is the x component of the sixteen control points of a surface patch. The 
matrix Q, 1s 


P, Ps. Po, Pis. 
Ro eee no, Sean 
S| Pewee peed aE 
P, Ps Pic Pie, 
and similarly for Q, and Q,. 

Since we must provide three 4x4 matrices, one for each of component x, y, 
and z, it can be seen that we have specified the 48 degrees of freedom as in the 
algebraic form. 

As can be seen by the above equations, a bicubic surface patch can be defined 
by a set of 16 control points and a basis matrix. By manipulating the control 


points, we change the shape of the surface as constrained by the basis matrix. 


We take this knowledge with us as we design our functions in the next chapter. 
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IV. DESIGN AND IMPLEMENTATION OF OUR SURFACE FUNCTIONS 


A. OVERVIEW 

Having the mathematics developed in the previous chapter does us no good if 
we can not to put it to use. In our case, this means being able to draw a 
parametric bicubic surface. It is at this point that we begin to see how the 
mathematics can be combined with the power of the computer and the graphics 
workstation. 

Up to this point. we have dealt with two forms of the parametric cubic curve 
and parametric bicubic surface - the algebraic form and the geometric form. The 
question now is which one shall we work with? 

Deciding what form to use depends largely on the application. If we are given 
or know the algebraic equations of the curve, then the reasonable choice is the 
algebraic form. If we plan to do surface fitting of data or interactive design the 
choice is the geometric. We choose to use the geometric form. Our primary 
reason for choosing it is that the geometric form offers us a greater insight into 
the control and behavior of curves and surfaces than is otherwise available with 
the classical algebraic formulation. It should be noted, however, that it is 
possible, through mathematical manipulation, to convert from one form to the 


other [Ref. 1l:pp. 164]. 
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B. METHODOLOGY 

How one proceeds to generate a surface impacts usefulness, flexibility, and the 
ability to understand. As we stated in the first chapter, the built-in functions 
provided by a graphics environment are not always exactly what we want. This 
is the case in the IRIS graphics environment. Although the IRIS provides support 
for bicubic surface patches. it does not support surface patch decomposition. 
What we want is the capability to do both. We also want this capability without 
sacrificing what a user already knows about how the IRIS supports bicubic surface 
patches. To achieve this. we have developed a set of parallel routines [Appendix 
A] that provide nearly all the functionality as the standard IRIS functions while 
at the same time providing the user with extended support via three additional 
functions. These extended functions allow the user to have access to the 
triangular polygons that our new functions generate during the construction of 
the bicubic surface patch. That is, the IRIS user is able to use the new functions 
in the same way as he would use the standard functions by substituting the 
names of the new routines in place of the standard IRIS routines. If, however, the 
user wishes to be able to have access to the individual triangular polygons that 


make up the surface, he has only three additional routines to learn. 


C. PARALLELING THE IRIS SUPPORTED FUNCTIONS 
The IRIS graphics environment has available five functions for defining and 


generating parametric bicubic surface patches. Those five functions are: 
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- defbasis defines a basis matrix 


- patchbasis sets the current basis matrices for both the s and the t parametric 
direction 


- patcheurves sets the number of curves used to represent a patch 
- patchprecision sets the precision at which the curves are drawn 


- patch draws the surface patch. 
A complete description of the functions and their arguments can be found in 
the IRIS Users Manual |Ref. 4]. 
To ease the pain of learning new functions, our new functions are 
syntactically identical to the standard IRIS functions with the exception that the 
new function names are the standard function names prefixed by the letter n. 


These parallel functions are: 


- ndefbasis 

- npatchbasis 

- npatcheurves 

- npatchprecision 
- npatch. 


The usage and the arguments of these parallel functions remain the same as 
the standard IRIS functions. The only difference that the user notices is that the 
wireframe drawn looks like a triangular mesh instead of the typical wireframe and 
that the function npatchprecision has no effect on the displayed image. 

While these routines seem to do what the old routines do, they are more 
powerful because they provide special extensions to the user. These extensions 


provide the capability to manipulate the surface patch as individual polygons. 
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D. TRIANGULAR DECOMPOSITION OF SURFACE PATCHES 
The extensions mentioned in the previous section are available to the user by 


using three additional routines: 


- Set User Routine for npatch provides an intercept function for handling the 
triangular polygons generated in the surface decomposition 


- User Routine provides the user a way to turn the intercept function on and 


off 


- Set Default Routine for npatch allows the user to return to the system 
defined intercept function. 


There is one argument to the function Set User Routine for npatch. This 
argument is the name of a user-defined function that expects to receive a 3x3 
array. This 3x3 array contains the three vertices of a triangle where each vertex is 
made up of an x, y, and z coordinate. The function User Routine expects one 
argument also. If this argument is zero, then the intercept function is turned off; 
i.e., the user's program cannot intercept the triangles composing the surface. 
Otherwise the function is activated allowing the user’s program to intercept the 
triangles comprising the surface. The function Set Default Routine for npatch 
does not expect any arguments. 

Using these functions, the user’s program has access to and can manipulate 
the individual triangular components of the surface patch. For example, an 
individual surface patch can be decomposed into triangular polygons and then via 
the user intercept function, each polygon can be subjected to an illumination 


model that produces a realistic looking surface in three-dimensional space. 
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The surface patch is decomposed into triangular polygons one at a time. For 
the users program to intercept these polygons, the program must have specified 
an intercept function via the Set User Routine for npateh and must have 
activated it via User Routine. Then, as each polygon is generated, the user's 
function can process them in any way desired. They can be stored, manipulated, 
altered. etc.. It is the user’s program that determines what to do with them. 
This feature provides a tremendous amount of flexibility, creativity, and 
applicability above what is currently available in the standard IRIS graphics 


environment support of surfaces. 


E. GENERAL GUIDELINES FOR USAGE 

To prevent any unnecessary problems in using our new functions, we need to 
establish a basic set of guidelines or sequences of events that should be followed. 
If the user does not want to use the special extensions, i.e. 
Set User Routine for npatch and User Routine, then a modified version of the 
standard IRIS setup steps for using surface patches can be followed. These steps 


are: 


define the appropriate curve bases using the ndefbasis function; 


- select a basis for the s and t parametric directions using the npatehbasis 
function; 


- select the number of curve segments to be drawn in each parametric direction 
using the npateheurves function; 


- draw the surface by using the npateh function. 


The only change to the standard IRIS setup is that it 1s not necessary to use 
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the npatchprecision function. This function does not effect the displayed image 
and its only purpose is to maintain consistency with the standard IRIS functions. 
If the user wishes to use the extensions via the Set User Routine for npatch 


and the User Routine functions, then the following steps must added: 


- An intercept function must be declared and defined before calling the 
Set User Routine for npatch function. This intercept function must be 
declared as a function returning an integer value (even though it is not used) 
and must be defined as receiving a 3x3 matrix of floating point numbers, 
where each row contains one set of x,y and z coordinates of an intercepted 
triangles vertex. The name of this intercept function will be the argument 
given to the Set User Routine for npatch function; 


- Activating/deactivating the intercept function via the User Routine function 
can be performed any time after the above step has been completed. 


By using these guidelines, a user should not have any difficulty in using these 
functions. As we will show in the next chapter, these functions are easy to use, 
efficient, and can provide some impressive results when a carefully chosen 


intercept function is used. 
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V. USAGE AND PERFORMANCE 


Having taken a brief tour of the functions we designed in the previous 
chapter. we need to provide some concrete examples of their usage, performance 


levels, and limitations. 


A. SAMPLE PROGRAMS 

Appendix B gives the listing for four sample programs using the new 
functions. Each program illustrates how the new functions can be integrated into 
the IRIS graphics programming environment. Program #1 draws 2 surface 
patches in wireframe representation. One surface is drawn using the standard 
IRIS functions while the other is drawn using the parallel functions. When this 
program is run the user notices the different appearance of the wireframe surface 
patch drawn with the new functions. It has the triangular mesh appearance 
described in the previous chapter. Program #2 shows how a _ user-defined 
intercept function can be used via the three extension functions we have designed. 
This program intercepts the triangles generated during the patches decomposition 
and puts them into an IRIS graphical object. Program #3 shows how the user- 
defined intercept function can be dynamically changed during execution and 
Program #4 shows how a well chosen intercept function can be used to produce 


realistic lighting of a surface. 
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B. PERFORMANCE COMPARISONS 

Knowing that our functions work, we would like to know how efficient they 
are. The way that we approach this question is to compare our new functions to 
the standard IRIS functions. Because we have designed our functions to be 
substitutable as a set for the standard functions, we do not have any problems in 
testing the relative performance of the sets. However, two words of caution are in 
order before comparing these two sets of functions. First, the IRIS implements 
many of its graphics primitives via special purpose hardware. This is the case 
with the function pateh. Therefore. its parallel function mnpateh, which is 
implemented in software is not as fast. Second. the new function npatchpreetsion 
does not affect the computation whereas the IRIS function patchprectsion does 
affect the computations. Keeping these points in mind, we have developed a 
simple benchmark program, listed in Appendix C, that we use to draw a 
wireframe representation of a surface patch 100 times. By executing this program 
10 times and getting the average times, we can get an idea of the performance of 
the parallel set of surface patch functions as compared to the standard IRIS 
surface patch functions. 

The way that we measure performance of a particular program is to use the 
UNIX time command. The time command returns, on program completion, the 
time in seconds for system time, user trme, and elapsed time. 

The benchmark program was executed in two different modes. In the first 


mode, the program was executed without the assistance of the IRIS’s floating 
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point accelerator (FFP) while in the second mode, the program was executed with 
the FFP. The results indicate that the standard IRIS functions are 350% faster 
without the FFP and 260% faster with the FFP. As expected, the new functions 
are slower, however considering that they were not designed to replace the IRIS 
functions. these results are good. Normally. one can expect an order of magnitude 


increase in performance when special hardware is used. 


C. LIMITATIONS 

Nothing that can be developed is without limitations. The reader should 
recall that it was certain limitations of the existing IRIS system that motivated 
this study. The functions we have designed and have implemented have allowed 
us to overcome certain limitations in the IRIS graphics environment. At the same 
time. these functions have their own limitations. 

The primary limitation of our parallel functions is speed. While these 
functions have been carefully implemented using efficient algorithms and data 
structures, they are not as fast as using special purpose hardware. Another 
limitation deals with the use of memory. The npatch function allocates memory 
to save each point on the surface patch. The number of points that are generated 
are proportional to the product of the desired number of curve segments in the s 
and t parametric directions. For example, to draw a surface patch with 10 curves 
in the s direction and 10 curves in the t directions requires at least enough 


memory to store 300 floating point numbers (one for each x, y, and z component). 
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To draw a 100x100 surface patch requires enough memory to store 30,000 floating 
point numbers. Assuming a floating point number requires 4 bytes, the 10x10 
patch requires 1.17 Kilobytes of memory, while the 100x100 patch requires 117.1 


Kilobytes of memory. 
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Reon eC TIONS POR FURTHER Ss PUDY 

Bicubic surface display and generation is an area of research in computer 
graphics that is exciting and important. Since the development of high- 
performance graphics systems, the demand for realism and real-time has increased 
significantly. Consequently. the need is great for continued creativity and 
exploration in the area. 

1. Development of Application Programs 

The power of these parallel surface functions we have created can only be 

derived through the use of the intercept functions. Whether they will be used to 
experiment with lighting and shading models or applied to fractal geometry can 
only be answered by time. However. it is through creative experimentation that 


these questions can be answered. Some areas for further work are: 


- Surface-fitting sampled data 

Surface-fitting is the process of constructing a representation to model the 
surface of an object based on a fairly large number of given data points. By 
taking these points and chosing an appropriate set of surface constraints, one 
can accurately reconstruct the surface. For example, during this work, we 
were given a set of digitized x, y and z coordinates for a human head. By 
successively extracting control points from the data, we were able to 
reproduce and display the head quite accurately. 
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- Data Reduction 

In many instances it is possible to reduce the amount of data needed to 
properly reconstruct a surface. For example, consider geographical terrain. 
Terrain that is relatively flat can be reconstructed with fewer surface patches 
that can mountainous terrain. The problem is that most terrain is sampled 
at discrete intervals, such as every 100 meters, whether it is flat or not. By 
applying some form of an Adaptive Subdivision Algorithm |Ref. 5] one can 
reduce the amount of primary and secondary storage while at the same time 
provide increased performance for display. 


- Lighting Models 
Because the user can intercept individual polygons comprising the surface, it 
is possible to subject each polygon to a lighting model. While we have 
provided a simple example of this, more sophisticated lighting models could 
be easily integrated through these parallel functions. 


- Realistic 3-D Objects 
The surfaces of many vehicles such as automobiles. aircraft, and ships can be 
constructed with bicubic surface patches. For example, constructing an 
object with surface patches and applying a lighting and shading model, one 
could develop an ship identification training system. Such a training system 
would be a valuable asset in military training environments, allowing the 
trainee to view a particular class of ship from any viewing angle. 
2. Improvement of Performance 
Real-time computer graphics requires efficient algorithms and data 
structures. While these functions were coded to be as efficient as possible, while 
preserving understandability, there is always room for improvement. One 
suggestion we have is to contact the developers of the IRIS graphics package for 
insights into improving our packages performance. Such contact may provide 


access to low-level graphic system routines and techniques that could dramatically 


improve performance. 
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Be CONGEUSIONS 

This study introduced the reader to the world of parametric bicubic surfaces. 
To do this. we provided some necessary definitions, terminology and mathematics. 
We also designed and implemented a set of software functions that take 
advantage of the information and given them to the reader for experimentation. 
The benefit that can be derived from the use of these functions can only be 


determined by the passage of time. 
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APPENDIX A —- FUNCTION SPECIFICATIONS 


NAME 
ndefbasis - defines a basis matrix 
SPECIFICATION 
ndefbasis(id, mat) 
long id: 
Matrix mat; 
DESCRIPTION 
ndefbasis allows the user to define basis matrices for use in the 
generation of patches. matriz is saved and is associated with itd. td 
may then be used in subsequent calls to npatchbasis. 
NAME 


npatchbasis - sets current basis matrices 


SPECIFICATION 


npatchbasis(sid, tid) 
long sid, tid; 


DESCRIP TION 
npatchbasis sets the current basis matrices (defined by ndefbasis } 


for both the s and t parametric directions of a surface patch. The 
current s and ¢ bases are used when the npatch command is issued. 
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NAME 
npatchcurves - sets number of curves used to represent a patch 
SPECIFICATION 


npatchcurves(scurves, tcurves) 
long scurves, tcurves; 


DESCRIPTION 


npatchcurves sets the current number of s and ¢ curves used to 
represent a patch as a wireframe. 


NAME 
npatchprecision - is a null function. 
SPECIFICATION 


npatchprecision(ssegments, tsegments) 
long ssegments, tsegments; 


DESCRIPTION 
npatchprecision has no functionality at the current time. It is used 


to maintain consistency with the standard IRIS function 
patchprecision. 
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NAME 
npatch - draws a surface patch 
SPECIFICATION 


npatch(geomx, geomy, geomz) 
Matrix geomx, geomy, geomz; 


DESCRIPTION 


npatch draws a surface patch using the current npatchbasis and 
npatchcurves. The shape of the patch is determined by the control 
points specified in geomz, geomy, and geomz. 


NAME 
Set User Routine for npatch - allows the user to specify an 
intercept function 
SPECIFICATION 
Set User Routine for npatch(fname) 
int (*fname) (); 
DESCRIPTION 


Set User Routine for npatch allows the user to set up a function 
that is capable of intercepting triangular polygons generated during 
the decomposition of a surface patch. The number of polygons 
generated is ( (scurves - 1) * (teurves- 1) * 2. ). 
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NAME 


User Routine - allows the user to turn the intercept function on and 


off 
SPECIFICATION 


User Routine(boolean) 
int boolean: 


DESCRIPTION 
User Routine acts like a switch allowing a user-defined intercept 
function to be turned on and off. By assigning boolean the value 0 the 


intercept function is turned off. Integer values other than O cause the 
intercept function to be turned on. 


NAME 


Set Default Routine for npatch - resets the intercept function 
to a system defined default 


SPECIFICATION 
Set Default Routine for npatch() 


DESCRIPTION 


Set Default Routine for npatch enables the user to choose the 
system defined intercept function poly(3, Triangle). 
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APPENDIX B — DEMONSTRATION PROGRAMS 


/* This is file "basis.h" */ 
#define HERMITE 0 
#define BEZIER 1 
#define CARDINAL 2 
#define BSPLINE 3 


/* the HERMITE BASIS MATRIX */ 
Matrix hermitematrix = { 

{20h e2 Osten On) 

{ -3.0, 3.0, ~2.0. -1.0 }, 

{ 0.0. 0.0. 1.00.0 }, 

{ 1.0. 0.0.0.0,0.0 } 


}: 


/* the CARDINAL BASIS MATRIX */ 
Matrix cardinalmatrix = { 

{ -0.5, 1.5, -1.5,0.5 }, 

{Ol 2-5 2-0450'5 }. 

{ -0.5, 0.0, 0.5.0.0}, 

{ 0.0, 1.0,0.0,0.0 } 


}; 


/* the BEZIER BASIS MATRIX */ 
Matrix beziermatrix = { 

{ -1.0, 3.0, -3.0. 1.0 }, 

{ 3.0, -6.0, 3.0,0.0 }, 

{ -3.0, 3.0, 0.0.0.0 }, 

{ 1.0.0.0,0.0,0.0 } 


}; 


/* the B-SPLINE BASIS MATRIX */ 
Matrix bsplinematrix = { 
{ -1.0/6.0, 3.0/6.0, -3.0/6.0, 1.0/6.0 }, 
{ 3.0/6.0, -6.0/6.0, 3.0/6.0, 0.0 }, 
{ -3.0/6.0, 0.0, 3.0/6.0,0.0  }, 
{ 1.0/6.0, 4.0/6.0, 1.0/6.0,0.0 } 
}; 
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/* This is file "geom.h" */ 


/* set up the geometry matrix of x coordinates * / 
Coord geomx[4][4] = { 

{ 0.0. 100.0, 200.0, 300.0 }. 

{ 0.0. 100.0, 200.0. 300.0 }. 

{ 1000.0. 900.0, 800.0, 700.0 }, 

{ 1000.0. 900.0, 800.0, 700.0 } 


re 


/* set up the geometry matrix of y coordinates */ 
Coord geomy[4][4] = { 

{ 400.0. 500.0, 600.0, 700.0 }, 

{ 0.0. 200.0. 400.0, 600.0 }, 

{ 0.0. 200.0. 400.0, 600.0 }, 

{ 400.0. 500.0, 600.0, 700.0 } 


1B 


/* set up the geometry matrix of z coordinates */ 
Coord geomz[4][4] = { 

{ 0.0. 200.0. 400.0, 800.0 be 

{ 0.0. 200.0. 400.0, 800.0 }, 

{ 0.0. 200.0, 400.0, 800.0 }, 

{ 0.0. 200.0, 400.0, 800.0 } 
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Program #1 


This program displays two wireframe images of the same 
surface patch. The patch drawn in the color YELLOW is 
produced by the standard IRIS patch functions and the 
patch drawn in the color RED is produced by the paralle] 
functions we have developed. 


One notices that the patch drawn via the parallel] functions 
has a triangular mesh appearance and that the call to the 
npatchprecision does not affect the displayed image as 

does the standard IRIS patchprecision function... 


Bee ee ee ee ee ee ee ee | ee ne aan 


#include "gi.h" /* IRIS graphics library */ 


#include "basis.h" 


#include “geom.h" 


#define S CURVES 10 
#define T CURVES 10 


main () 


{ 


/* Loop variables */ 
int pl, p2; 


/* initialize the graphics system */ 
ginit(); 

doublebuffer();: 

gceonfig(); 

cursoff(); 


/* set up the viewing parameters * / 
ortho(0.0, 1023.0, 0.0, 1023.0, -1023.0, 1023.0); 
viewport(0, 1023, 0, 767); 


/™ clear the graphics screen to BLACK */ 
color(BLACK); 


clear(); 


/" Associate an id number with a basis matrix * / 
defbasis(BEZIER, beziermatrix); 
defbasis(CARDINAL, cardinalmatrix); 
defbasis(BSPLINE, bsplinematrix); 


ndefbasis(BEZIER, beziermatrix); 


ndefbasis(CARDINAL, cardinalmatrix); 
ndefbasis(BSPLINE, bsplinematrix); 
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/** Specify how many curves in each 
parametric direction. **/ 
patchcurves(S CURVES .T CURVES): 


npatchcurves(S CURVES .T CURVES); 


/** Make the basis matrices different 
for each parametric direction. **/ 

patchbasis(BEZIER, CARDINAL); 

npatchbasis(BEZIER, CARDINAL); 


/** Cycle through the patch changing the 
precision that the individual curves 
comprising the patch are drawn. **/ 
for(pl = 10, p2 = 100; p1 < 100; pl += 5, p2 -= 5) { 


/* Draw the image via the IRIS functions */ 
viewport(0. 511, 0. 767): 

color(BLACK): 

clear(): 

color(YELLOW);: 

patchprecision(p1, p2): 

patch(geomx, geomy, geomz): 


/* Draw the image via the parallel functions */ 
viewport (512, 1023, 0. 767); 

color(BLACK); 

clear(); 

color(RED); 

npatchprecision(pl. p2); 

npatch(geomx, geomy, geomz); 


/* display the wireframe images */ 
swapbuffers(): 
sleep(1): 


/* clear the screen */ 
color(BLACK); 
clear(); 


} 


/* clear the graphics screen and exit */ 
color(BLACK); 


clear(); 
gexit (); 
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Program #2 


This program illustrates how the triangles formed during 
a surface patch decomposition can be put into an IRIS 
graphical object and subsequently displayed. 


FO ee ee ee ee ee en eet eee eae aaeeee nae 


#include "ol.h" /* IRIS graphics library = 


f#finclude "basis.h" 
#include "geom.h" 


#define ON 1 
#define S CURVES 10 
#define T CURVES 10 


/* Where we put our intercepted triangles */ 
Object intercepted object: 


main () 

{ 
/* declare the intercept function */ 
int intercept function(); 


/* initialize the graphics system */ 
ginit(); 

doublebu ffer(); 

gconfig(): 

cursoff(); 


/* Make the intital object. */ 
makeobj(intercepted object = genobj()); 
closeobj(); 


/* set up the viewing parameters * / 
ortho(0.0, 1023.0, 0.0, 1023.0, -1023.0, 1023.0); 
viewport (0, 1023, 0, 767); 


/* clear the graphics screen to BLACK */ 
color(BLACK); 
clear(); 


/* Associate an id number with a basis matrix * / 
ndefbasis(BEZIER, beziermatrix); 
ndefbasis(CARDINAL, cardinalmatrix); 
ndefbasis(BSPLINE, bsplinematrix); 


/** Specify how many curves in each 
parametric direction. ** / 


npatchcurves(S CURVES , T CURVES): 
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/** Make the basis matrices different 
for each parametric direction.** / 


npatchbasis(BEZIER, BSPLINE); 


/** Set up an intercept function to grab 
the triangles from the surface patch 
and turn it on. **/ 
Set User Routine for npatch(intercept function); 


User Routine(ON): 


/* Call the npatch function. */ 
npatch(geomx. geomy, geomz); 


/* Display the surface patch. */ 
color(YELLOW): 
callobj(intercepted object); 
swapbuffers(): 

sleep(10); 


/* clear the graphics screen and exit */ 
color(BLACK): 

clear(): 

gexit(); 


Me ak whe abe abe ade abe ale ade afc ade abe ade abe ale ade abe abe ake ah abe ale ale ae ade ae ae ade ade ae ade ahead ae ae ade ade ae ae ae abe ade abe eae ae abe eae hehe eee oe 


This is the user-defined intercept function that 
handles the individual triangles generated 
during the decomposition of a surface patch. 


eee ee eee renee ee ee eM / 


int intercept function(triangle) 

float triangle[3][3]; 

{ 
/* Open up the object and put in the triangles */ 
editobj(intercepted object): 
move (triangle[0]/0], triangle(0][1], triangle{0] [2] 
draw(triangle(1][0], triangle/1][1], triangle[1][2] 
draw(triangle(2]/0], triangle([2][1], triangle[2][2] 
draw/(triangle(0][0], triangle[0][1], triangle[0][2] 
closeobj(); 


3 


3 


Sd 


$ 
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Program #3 


This programs illustrates how the intercept function 
can be changed as the program runs. 


we We ee ee OO ee ee ees 


f#include "gl.h" /* IRIS graphics library */ 


#include "basis.h" 
#include "geom.h" 


/** Define a new type that is a pointer to a function 
that returns an integer. **/ 
typedef int (“Function pointer) () : 


#define ON 1 
#define S CURVES 10 
#define T CURVES 10 


main() 

{ 
/* Loop variable */ 
int count; 


/* Declare an array of pointers to functions */ 
Function pointer intercept functions/4]; 


/* declare the intercept functions */ 
int intercept function1(): 
int intercept function2():; 
int intercept function3(); 
int intercept function4() 


§ 


/** Initialize the array of pointers to the 
intercept functions. ** / 

intercept functions[0] = intercept function]; 

intercept functions({1] = intercept function?2; 

intercept functions[2] = intercept function3; 

intercept functions/3] = intercept function4; 


/* initialize the graphics system * / 
ginit(); 
doublebuffer(); 


gconfig(); 
cursoff(); 


/* set up the viewing parameters * / 


ortho(0.0, 1023.0, 0.0, 1023.0, -1023.0, 1023.0); 
viewport (0, 1023, 0, 767); 
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/* clear the graphics screen to BLACK */ 
color(BLACK); 
clear(); 


/* Associate an id number with a basis matrix */ 
ndefbasis(BEZIER, beziermatrix); 
ndefbasis(CARDINAL. cardinalmatrix); 
ndefbasis(BSPLINE. bsplinematrix); 


/** Specify how many curves in each 
parametric direction. **/ 


npatchcurves(S CURVES .T CURVES); 


/** Make the basis matrices different 
for each parametric direction. ** / 


npatchbasis(BSPLINE, CARDINAL); 


/* Initially use the default intercept function */ 
Set Default Routine for npatch(): 
User Routine(ON): 


/* Step through each intercept function * / 
for(count = 0; count < 4: count+-+) { 


/** Set up an intercept function to grab 
the triangles from the surface patch. **/ 
Set User Routine for npatch(intercept functions[count)); 


/** Call the npatch function using the current 
intercept function. **/ 
npatch(geomx, geomy, geomz): 


/* Display what the intercept function did. */ 
swapbuffers(); 
sleep(2); 


/* Clear the screen and do another one. */ 
color(BLACK); 
clear(); 
} 
/* clear the graphics screen and exit */ 
color(BLACK); 
clear(); 
gexit(); 
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This intercept function draws each triangle RED. 


+ 0S Ee OO Oe Oe ee ee eee | 


int intercept function1(triangle) 
float triangle|3]/[3]; 
{ 

color(RED): 

poly(3, triangle): 


De eee ee eee ee ee ee eee 


This intercept function draws each triangle YELLOW. 


Fee ee ee ee ee eee ee Ye armada et edie eset cs aaa | 
int intercept function2(triangle) 
float triangle/3][3]; 


{ 
color(YELLOW): 


poly (3, triangle): 


| haclatiadiaciaciaciatiatiaciatiadiadiacn taunts ndiatadateiiadan i 


This intercept function draws each triangle GREEN. 


Cn ee tg ie ii iaiee ae eeeeneedis coe ae 
int intercept function3(triangle) 
float triangle([3}[3]; 


{ 
color(GREEN): 


poly(3, triangle): 


Pe ee eee ee SS ton oe a eae ec ee 


This intercept function draws each triangle BLUE. 
{Fee eA eOeA eC REE OEE HEE SEER ERE ERED ERE O DLE OO EEE OH / 


int intercept function4(triangle) 
float triangle/[3][3]; 


{ 
color(BLUE); 


poly(3, triangle); 
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Program #4 


This program illustrates how the intercept function 
can be combined with an illumination model to 
provide a realistic, illuminated, three-dimensional 
solid-filled curved surface. 


Se eT ee ee oe” ee a / 


#include "gl.h" /* IRIS graphics library */ 


#include "device.h" 
#include "basis.h" 


#include "geom.h" 


#define ON 1 
#define S CURVES 25 
#define T CURVES 25 


main () 


{ 


/* Declare the intercept function. */ 
int light poly(): 


/* loop variables * / 
int i,j: 


/* Initialize the graphics system. */ 
ginit(): 

singlebuffer (): 

gconfig(); 

cursoff(); 


/* Clear the display. */ 
color(BLACK); 
clear(); 


/* Set up new viewing parameters. */ 
ortho(0.0, 1023.0. 0.0, 1023.0, -1023.0, 1023.0); 
viewport(0, 1023, 0, 767); 


/* Clear drawing area. */ 
color(CYAN); 
clear(); 


/* Set for hidden surface elimination. * / 
setdepth(Ox3F FF, 0xC000); 

zclear(); 

zbuffer(TRUE); 
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/* Load the color map ramp with a grey scale. */ 
for (i= 0;i< 256; i++) 


} 


mapcolor(8-+i, i, i, i); 


/* Associate an id with a basis matrix. */ 
ndefbasis(BEZIER, beziermatrix): 
ndefbasis(CARDINAL, cardinalmatrix); 
ndefbasis(BSPLINE, bsplinematrix); 


/** Provide a different basis for each 


parametric direction. **/ 


npatchbasis(CARDINAL. BEZIER); 


/** Provide the number of curves in 


each parametric direction. **/ 


npatchcurves(S CURVES .T CURVES): 


/* Set up the intercept function. */ 
Set User Routine for npatch(light poly): 
User Routine(ON); 


while(TRUE) 


{ 


/* Clear the z buffer. */ 
zclear(); 


/* Hold display if MOUSE2 is down. */ 
if(getbutton(MOUSE2)) 


/* Resume when MOUSEL is pressed */ 
while(!getbutton(MOUSE1)) ; 


} 


/* Exit when MOUSE1 and MOUSE2 and MOUSES3 are down. */ 
if(getbutton(MOUSE1) && getbutton(MOUSE2) && getbutton(MOUSE3)) 
break; 


/* Clear the drawing area. */ 
color(CYAN); 


clear(); 


/* Set the current color. */ 
color(BLACK); 


/* Draw initial surface patch. */ 
npatch(geomx, geomy, geomz); 
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/** Change the y coordinates to get a 
different surface patch the next 
time we draw display. **/ 

for (i=1; 1<3; i++) 

{ 


for (j=0: j<4; j++) 
geomy/(i][{j] = (float) (lrand48() % 900); 


} 


/* Draw the surface patch. */ 
npatch(geomx, geomy, geomz): 


pun Change the y coordinate values to 
make another surface patch. ** / 

for (i=1; 1<3; i++) 

{ 


for (j=0: j<4; j++) 


geomy/|i][j] = (lrand48()%2)*geomy(i][j]: 


} 


/* Clean up and exit the program. */ 
color(BLACK); 

clear(): 

gexit (): 
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The user-defined intercept function used to grab the 
triangles generated in the surface patch decomposition. 
Se i ee eee RS YL 
int light poly(Triangle) 
Coord Triangle|3]{3]; 
{ 
/* Put each triangle through an illumination model. */ 
lightthepoly (Triangle, 3, 350.0, -1750.0, 350.0, 350.0, 1750.0, 350.0, 9, 264); 


ol 
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hghtpoly.c 


It is a routine that computes lighting for a polygon based 
upon the angle between the Normal vector of the polygon 
and the direction to the light source. 


light thepoly(xyz,ncoords,ax,ay,aZ,]x,ly.1z,colormin,colormax) 
xyz{][3] = floating coords of the polygon. 
ncoords = number of coordinates. 


ax,ay.az = interior point of the whole object. Used to determine 
outward facing normal of the polygon. This is the same 
point of reference that would be used for backface 
polygon removal. 


Ix,ly.lz = vector pointing im direction of the lhght source. 


colormin, colormax = indices used for the colors assigned to this 
polygon. The user is responsible for setting 
up the color ramp. 


Note: the routine also puts the polygons out ordered counterclockwise 
with respect to the interior point for ease of backface polygon 
removal. 


A I A ee end cee re eee ee ell 


#include <gl.h> 
#imnclude <math.h> 


#define PIDIV2 1.570796327 
#define CLOCKWISE 1 
#define ROW 3 


lightthepoly(xyz,ncoords,ax,ay,az,1x,ly,1z,colormin.colormax) 
Coord xyz(][3]; 
unsigned int ncoords; 


Coord ax,ay,az: /* interior point of the whole object. */ 
Coord Ix,ly,lz;_ /* direction to the light source sa 
int colormin,colormax; /* color min/max indices “*/ 


{ 


/* temp coord hold */ 
Coord *txyz; 


/* loop temps */ 
register unsigned short int i,j; 
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/* direction test function */ 
int npoly orient(); 


/* vectors used to compute the polygon’s normal */ 
Coord v1(3]}.v2[3]: 


/* the polygon’s normal * / 
Coord normal|3}: 


/* normal’s magnitude */ 
Coord normalmag; 


/* light’s magnitude */ 
Coord lightmag; 


/* dot product of N and L */ 
double dotprod: 


/* angle between N and L */ 
float radians; 


/* color to use in drawing the polygon */ 
unsigned short int colortouse; 


/* allocate memory for a temporary array */ 
txyz = (Coord *) calloc((ncoords * 3), sizeof(Coord) ): 


/** orient the polygon so that its counterclockwise with respect 
to the interior point **/ 
if(npoly orient (ncoords,xyz,ax,ay,az) == CLOCKWISE) 


/* the polygon is clockwise. reverse it. * / 
for(i=0: i < ncoords; i=i+1) 
{ 

for(j=0; j < ROW; j=j+1) 

{ 


} 
} 
} 


else 


*(txyz + (i* ROW) + j) = xyz|ncoords-i-1]|j]; 


/* no need to reverse */ 
for(i=0; i < ncoords; i=i+1) 
{ 
for(j=0; j < ROW; j=j+1) 
{ 


} 
} 
} 


*(txyz + (ROW * i) + j) = xyz[i][i); 
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/* the coordinates are ordered counterclockwise in array txyz */ 


/** compute the normal vector for the polygon using the first 
three vertices... **/ 


/* compute the first vector to use in the computation */ 

v1[0] = *(txyz + 6) - *(txyz + 3); /* txyz[2][0] - txyz1}[0] */ 
vi[1] = *(txyz + 7) - *(txyz + 4); /* txyz[2][1] - txyz(t][1] */ 
vi[2] = *(txyz + 8) - *(txyz + 5); /* txyz[2][2] - txyz[1][2] */ 


/* compute the second vector to use in computing the normal */ 
v2[0] = *(txyz ) - *(txyz + 3); /* txyz[0][O] - txyz[1][0] */ 
v2[1] = *(txyz + 1) - *(txyz + 4); /* txyz[0][1] - txyz[4][1] */ 
v2[2] = *(txyz + 2) - *(txyz + 5); /* txyz[0][2] - txyz[1][2] */ 


/* the normal is v1 x v2 */ 

normal(0] = v1[1]*v2[2] - v1[2]*v2/1]: 
normal(1] = v1[2]*v2(0] - v1[0|*v2[2]: 
normal({2} = v1[0]*v2([1] - v1[1]*v2/[0}]: 


/* compute the magnitude of the normal */ 
normalmag = sqrt((normal(0|*normal[0]) + 
(normal({1]*normal{1])+ 
(normal |2]*normal({2})); 


/* compute the magnitude of the light * / 
lightmag = sqrt( (Ix * Ix) + (ly * ly) + (Iz * Iz)); 


/* compute N .L (normal dot product with the light source direction) */ 
dotprod = (double) ((normal{0] * 1x) 

+ (normal|1] * ly) 

+ (normal{2] * Iz)); 


/* compute the unit normal */ 
dotprod = (double) ((dotprod/(normalmag * lightmag))); 


/* dotprod = cos(theta) of the angle between N and L. 


Convert to angle in radians */ 
radians = acos(dotprod); 
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/* compute the color we should use */ 
if(-PIDIV2 <= radians && radians <= PIDIV2) 


{ 


/* if the angle is negative, set to positive */ 
if(radians < 0.0) 


radians = -radians; 


} 


colortouse = ({colormax-colormin) /PIDIV2)*(PIDIV2-radians)+colormin: 


} 


else 


{ 


colortouse = colormin: 


} 


/* set the color */ 
color(colortouse); 


/* draw the poly */ 
polf(ncoords.txyz): 


/* free up memory allocated for the temporary array */ 
cfree(txyz, (ncoords * 3), sizeof(Coord)); 
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npoly orient.c 


This routine determines a polygon’s orientation 
with respect to its normal and a reference point. 
Orientation is either clockwise or counter-clockwise. 
The point of reference must not lie in the polygon’s 
plane. 


aalaciaiaciacia inci ite aidan aciadilad oda das aad cd sc lag aces Sy) NI The 


#include <gl.h> 
f#include <math.h> 


int npoly orient(ncoords.xyz,xinside,yinside.zinside) 
unsigned int ncoords; 
Coord xyz(][3]; 
Coord sxonside, yinside, zinside: 
{ 
/* loop temps */ 
register unsigned short int i,j: 


/* center coordinate of the polygon */ 
Coord center{3}; 


/** vector hold locations for the vectors that run 
from the center coordinate to the points of the 
polygon ** / 

Coord a[3], b[3]; 


/** points on line containing normal that are 
on opposite sides of the plane containing 
the polygon. **/ 

Coord xn{3], xmn/{3}; 


/* distance to point n from pt inside. */ 
float distton; 


/* distance to point -n from pt inside. */ 
float disttomn; 


/* the normal vector computed from a x b */ 
Coord normal([3}; 


/* compute the center coordinate of the polygon */ 
center/0] = 0.0; 
center/1] = 0.0; 
center[2| = 0.0; 


96 


for (i=O0; 1 < ncoords; i++) 


{ 


for(j=0: j < 3; j++) 
center|j] += xyz[i][}j]; 
} 


/* divide out by the number of coordinates */ 
for (j=0; j < 3; j++) 


center(j] = center(j] /(float)ncoords; 


/** check the first 2 coordinates of the 
polygon for their direction **/ 


/** compute vector a. It runs from the 
center coordinate to coordinate 0 **/ 
for(j=0; j < 3: j++) 


aij] = xyz(0][j] - center|j]; 


/** compute vector b. It runs from the 
center coordinate to coordinate 1 **/ 


for (j=0; j <3; j++) 


b[j] = xy2[1][j] - center[j]; 


/* compute a x b to get the normal vector */ 
normal[0] = a[1|*b[2] - a(2]*b[1]; 
normal|1] = a[2]*b(0] - a[0|*b/[2]; 
normal(2]| = a[0|*b[1] - a[1]*b[0]; 


/** compute point n, offset pt from center in 
direction of normal **/ 


for (j=0; j < 3; j++) 


xn(j] = center(j] + normal(j]; 


/** compute point -n, offset pt from center 
in opposite direction from normal. ** / 

for (j=0; j < 3; j++) 

{ 


xmn|[j] = center|j] - normal|j]; 
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/* compute the distance the inside pt is from point n */ 
distton = sqrt((xn[O] - xinside) * (xn[0] - xinside) + 
(xm{1] - yinside) * (xm[1] - yinside) + 
(xn{2] - zinside) * (xn[2] - zinside)); 


/* compute the distance the inside pt is from point -n */ 

disttomn = sqrt((sxann/[0] - xinside) * (xamn|[0] - xinside) + 
(sxa¢mn[1] - yinside) * (xann[1] - yinside) + 
(xamn/[2] - zinside) * (xann[2] - zinside)); 


/** if the dist(n) < dist(-n). then n points back towards the 
inside point and is on the same side of the plane as inside. 
a x bis then clockwise. **/ 

if(distton < disttomn) 


return(1); /* clockwise */ 


} 


else 


{ 


return(0); /* counterclockwise */ 


} 
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APPENDIX C — BENCHMARK PROGRAM 
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This is the BENCHMARK PROGRAM used to test the 
relative performance of the standard IRIS functions 
and the parallel routines to those standard 

functions. 


me me mi ok mie sie rie aie oe rie i ahs aie ie ake sic ica aici eae aie ic aie ioe he ake eae ie ica akc oe ie aie aie i eae i ic aie aie ac eee aie sic ae oe aie oe 7 oh 


#include "gl.h" /* IRIS graphics library */ 


#include "basis.h" 
#include "geom.h" 


#define 


#define 
# define 
#define 
#define 


main() 


{ 


IRIS /* Which set to test switch. */ 


MAX TIMES THRU 100 
S CURVES ~ 25 
T CURVES ~ 25 


ONE 1 


int times thru: 


/* initialize the graphics system */ 
ginit(): 
doublebuffer(): 


gconfig(): 
cursoff(); 


/* clear the graphics screen */ 
color(BLACK); 
clear(); 


/* set up the viewing parameters * / 
ortho(0.0, 1023.0, 0.0, 1023.0, -1023.0, 1023.0); 
viewport(0, 1023, 0, 767); 


/* clear the graphics screen to CYAN */ 
color(CYAN); 


clear(); 
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#ifdef IRIS 
/* Use the standard IRIS functions */ 
defbasis(BEZIER, beziermatrix); 
defbasis(CARDINAL. cardinalmatrix); 
defbasis(BSPLINE, bsplinematrix); 


patchbasis(BEZIER, BEZIER); 
patchcurves(S CURVES .T CURVES); 
patchprecision(ONE, ONE); 

f#else 
/* Use the parallel functions */ 
ndefbasis(BEZIER, beziermatrix); 
ndefbasis(CARDINAL, cardinalmatrix); 
ndefbasis(BSPLINE, bsplinematrix); 


npatchbasis(BEZIER, BEZIER); 

npatchcurves(S CURVES .T CURVES): 

npatchprecision(ONE, ONE): 
#endif 
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for(times thru = 0: times thru < MAX TIMES THRU: times thru+-+) { 


/* clear the graphics screen to CYAN each time thru */ 
color(CYAN); 
clear(); 


/* draw the wireframe surface patch in BLACK */ 


color(BLACK): 
#ifdef IRIS 
/* Using the IRIS function */ 
patch(geomx, geomy, geomz); 
#else 
/* Using the parallel function */ 
npatch(geomx, geomy, geomz); 
#endif 


} 


/* display the wireframe image */ 
swapbuffers(); 


/* clear the graphics screen and exit */ 
color(BLACK); 

clear(); 

gexit(); 
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APPENDIX D — PARALLEL FUNCTIONS SOURCE CODE 


ee eee ee ee akg i ee ee eae 


FILE: npatch.c 

Author ..... Gary W. TAYLOR (Captain USMC) 
Date 2. 1 December 1986 

Place®...-. Naval Postgraduate School, Monterey CA 
Environment ..... 


Silicon Graphics, Inc., IRIS 2400 
graphics workstation, UNIX operating 
system (GL2-W3.4). 


The following C source code provides a set 
of "shadow" routines to parallel the standard 
IRIS 2400 graphics workstation surface patch 
routines. These parallel routines provide 
their user the capability to generate 
solid-filled parametric bicubic surface patches. 


As of the current date, there are no known 
side-effects or bugs associated with using 
these functions. 


Limitations ..... 


These functions can be used in immediate mode only. 


ee eee ee ee ee ee alee ee cede ee eee 
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#include "gl.h" /* IRIS graphics library */ 
#include "stdio.h" 


#define OFF 0 


/* Here is how we define a surface point. */ 
typedef struct { 

Coord x: 

Coord y: 

Coord z: 
} Point: 


/* Structure used to track user supplied basis matrices. */ 
static struct list elem { 


/* Matrix id number. wy 
long nid: 


/* Pointer to the basis matrix. */ 
float *nmatrix; 


/* Pointer to the next basis matrix. */ 
struct list elem *next elem ptr; 


he 
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/* Used in forward difference computation in the U direction. */ 
static Matrix Precision Matrix U; 


/* Used in forward difference computation in the V direction. */ 
static Matrix Precision Matrix V; 


/* The default intercept function. */ 
static int Default system routine (): 


/* Pointer to the currently defined U basis matrix. */ 
static oat *current U_ basis: 


/* Pointer to the currently defined V basis matrix. */ 
static float *current V_ basis: 


/* How many curves in the U direction. */ 
static long UCURVES:; 


/* How many curves in the V direction. */ 
static long VCURVES; 


/* How may curve segments in the U direction. */ 


static long USEGMENTS; 


/* How may curve segments in the V direction. */ 


static long VSEGMENTS; 


/* Pointer to linked list of user supplied basis matrices.* / 
static struct list elem *head of list = NULL; 


/* Set initial user routine to the default. */ 
static int (* User routine) () = Default system routine; 


/* Initially do not call user’s function. */ 
staticint User routine is on = OFF; 
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_Default system routine 


This routine is what the system automatically does if the 
user does not supply a particular intercept function. 
ns ea a a i te eee Ree ee 


static int Default system routine (Triangle) 
Coord Triangle[3][3]; 


{ 


/* Draw the polygon with a standard IRIS function. */ 
poly (3, Triangle); 


ee ee NES RR ERE EEE ee ee 


Set User Routine for npatch 


This function allows the user to supply an intercept function 
to be used to manipulate the triangles that are generated in 
the decomposition of a surface patch. 


i Ce ee ee ee eS 


void Set User Routine for npatch (routine) 
int (*routine) (); 


{ 


/* Save the pointer to the user’s function. */ 
_User routine = routine; 
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User Routine 


This function allows the user to turn the intercept function 
on or off at will. 


ae ee eece ne ne ee eg nn 


void User Routine (boolean) 
int boolean; 


Yael 
If boolean = 0 then the intercept function is turned OFF. 
If boolean != 0 then the intercept function is turned ON. 


cf 


User routine is on = boolean: 


Dee ee eee eee ee ee eee ee eee ene 


Set Default Routine for npatch 


This function allows the user to reset the intercept routine 
to the same routine used by the system. 


ee ee ee ee ee oe memiee — ) ciimiemmiermei gia eke 


void Set Default Routine for npatch () 


{ 


/* Point to the default intercept function. */ 
_User routine = Default system routine; 
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ndefbasis 


This function is equivalent to the IRIS defbasis function 
in that it allows the user to define a basis matrix for 

use in the generation of patches. matrix is saved and 

is associated with id. id may then be used in subsequent 


calls to npatchbasis. 
ihe t ee ee ie | 
void ndefbasis (id, matrix) 


long id: 
Matrix matrix: 


{ 


/* Special processing first time this function is called. */ 
static int has been called = FALSE; 


/* Data structure pointer for new entry. */ 
struct list elem *new elem to add: 


/* Pointer to search linked list. * / 
struct list elem *walking ptr: 


/* Pointer to a copy of the user supplied basis matrix. * / 


float *pmatrix: 


/* Loop variables */ 
int row: 
int column: 


/* Get memory for the new data elements. */ 
new elem to add = (struct list elem *) malloc (sizeof (struct list_elem)); 
pmatrix = (float *) calloc (sizeof (Matrix), sizeof (float)); 


/* Make a copy of the basis matrix passed in by the user. */ 


for (row = 0: row < 4; row+-+) 
for (column = 0; column < 4; column+-+) { 


*(pmatrix + (4 * row) + column) = matrix{row]/column]; 


} 


/* Associate the user supplied id to this basis matrix. */ 
new elem to add -> nid = id; 


/* Point to the copied basis matrix. * / 
new elem to add -> nmatrix = pmatrix; 
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/** Determine how to add this information into the linked list 
of basis matrices. **/ 


/* Does a list already exist? */ 
switch (has been called) { 


case TRUE: 


/* Point to the beginning of the list. */ 
walking ptr = head of list; 


/** Traverse the list looking to see if the id number 
already exists. **/ 
while ((walking ptr -> nid != id) && 
(walking ptr-> next elem ptr != head of list)) { 


/* Walk through the linked list. */ 
walking ptr = walking ptr-> next elem ptr; 


} 


/* id already exists so we can reuse its allocated memory. */ 
if (walking ptr -> nid == id) { 


/* Get rid of the old basis matrix. */ 
cfree (walking ptr -> nmatrix, sizeof (Matrix), sizeof (float) ); 


/* Point to the replacement matrix. */ 
walking ptr -> nmatrix = pmatrix; 


/* Get rid of the un-needed data structure. */ 
cfree (new elem to add, 1, sizeof (struct list elem)); 
else { /* The id does not exist */ 
/** Manipulate the pointers to add the new 
data element to the linked list. **/ 
new elem to add -> next elem ptr = head of list -> next elem ptr: 


head of list > next elem ptr = new elem to add; 


} 


break; 
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case FALSE: 
/* No linked list of basis matrices exists so we start up one. */ 
/* Create the pointer to the front of the list. */ 

head of list = new elem to add; 


head of list -> next elem ptr = head of list: 


/* Make sure we do not do this again. */ 
has been called = TRUE; 


break; 


default: 


break; 
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npatchbasis 


This function is equivalent to the IRIS patchbasis function 


in that it sets the current basis matrices for both the 
U and V parametric directions of a surface patch. The current 


U and V bases are used when the npatch command is issued. 


Tee ee ee ee ee ee ee ee ee eS ee ee 


int npatchbasis (uid, vid) 
long uid, vid; 


{ 


struct list elem *walking ptr; 


/* ERROR: no linked list of basis matrices. */ 
if (head of list == NULL) { 


fprintf (stderr, "Opatchbasis: no basis matrices definedO); 


exit (-1); 


} 
/* Traverse the list looking for the desired U basis matrix. */ 


walking ptr = head of hist; 
while ((walking ptr -> nid != uid) && 
(walking ptr -> next elem ptr != head of list)) { 


walking ptr = walking ptr -> next elem ptr; 


} 


if ((walking ptr -> nid != uid) && 
(walking ptr -> next elem ptr == head of list)) { 


/* ERROR: U basis matrix does not exist in the linked list. */ 
fprintf (stderr, "Opatchbasis: undefined U basis matrix %d0, uid); 
exit (-1); 


} 


current U_ basis = walking ptr -> nmatrix; 
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/* Traverse the list looking for the desired V basis matrix. */ 
walking ptr = head of list: 
while ((walking ptr -> nid != vid) && 

(walking ptr -> next elem ptr != head of list)) { 


walking ptr = walking ptr-> next elem ptr; 


} 


if ((walking ptr -> nid != vid) && 
(walking ptr -> next elem ptr == head of list)) { 


/* ERROR: V basis matrix does not exist in the linked list. */ 


fprintf (stderr, "Opatchbasis: undefined V basis matrix %dO. vid); 
exit (-1): 


} 


current V_ basis = walking ptr -> nmatrix; 
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npatchcurves 


This function is similar to the IRIS patchcurves command. 
ucurves and vcurves set the subdivision parameters used in 
decomposing the surface patch. An individual patch will be 
decomposed into a (ucurves -1) * (veurves -1) grid with 
each grid generating two triangular polygons. 


Fe ee ee ee ae eset en ent oly, rear ae a 


void npatchcurves (ucurves, vcurves) 
long ucurves, vcurves; 


{ 


/* Prevent an inappropriate number of curves. */ 
UCURVES = (ucurves < 2 ? 2: ucurves - 1); 
VCURVES = (vcurves < 2 ? 2: veurves - 1); 


/** Set up the Precision Matrix U used in the forward difference 
along the U direction. **/ 


Precision Matrix U[0][0] = 6.0 / (float) (UCURVES * UCURVES * UCURVES); 
Precision Matrix _U[1][0] = 6.0 / (float) (UCURVES * UCURVES * UCURVES); 


Precision Matrix U(1][1] = 2.0 / (float) (UCURVES * UCURVES); 


Precision Matrix U[2][0] = 1.0 / (float) (UCURVES * UCURVES * UCURVES): 


Precision Matrix U[2][1] = 1.0 / (float) (UCURVES * UCURVES); 
Precision Matrix U[(2][2] = 1.0 / (float) (UCURVES): 
Precision Matrix U(3][3] = 1.0; 


/** Set up the Precision Matrix V used in the forward difference 
along the V direction. **/ 


Precision Matrix V(0][0] = 6.0 / (float) (VCURVES * VCURVES * VCURVES); 
Precision Matrix V[1][0] = 6.0 / (float) (VCURVES * VCURVES * VCURVES): 
Precision Matrix V[1][1] = 2.0 / (float) (VCURVES * VCURVES); 

Precision Matrix V(2][0] = 1.0 / (float) (VCURVES * VCURVES * VCURVES); 
Precision Matrix V(2][1] = 1.0 / (float) (VCURVES * VCURVES); 

Precision Matrix V([2][2] = 1.0 / (float) (VCURVES); 

Precision Matrix V(3][3] = 1.0; 
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npatchprecision 


This function is similar to the IRIS patchprecision routine 
but is only used to maintain consistency with the IRIS 
routines. Its results are not used by any other function 

in the suite. 


ieee en ee ee ee ee ee eee 


void npatchprecision (usegments, vsegments) 
long usegments. vsegments: 


{ 


/* Prevent an inappropriate number of segments. */ 
USEGMENTS = (usegments < 2 ? 2 : usegments - 1); 
VSEGMENTS = (vsegments < 2 ? 2: vsegments - 1); 
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nspeckle 


This function finishes computation and hands off each triangle 
to an intercept function. 


aa i achat cada is ve nada a dag sg od 


static void nspeckle (Coord array) 
Point * Coord array(4]; 


{ 


register Point * Patch array; 


/* Get enough memory to hold all the points that will be generated. */ 
Patch array = (Point *) calloc ((UCURVES + 1) * (VCURVES + 1), 
sizeof (Point));: 


{ 


/* Pointer to a particular point. */ 
register Point * Where; 


register unsigned int total points; 
register unsigned int point count; 
register unsigned int t count; 
register unsigned int count; 
register unsigned int Row; 
register unsigned int Column; 


/* Used in generating points on the surface. */ 
Matrix control matrix; 


/* Intermediate matrix to hold mathematical results. */ 
Matrix interI; 


/* For every point in the the U parametric direction. * / 
for (point count = 0, total points = 0; 


point count <= UCURVES; point count+-+) { 


/* Build a control matrix for the current curve. */ 
for (count = 0; count < 4; count+-+) { 


Where = (Point *) (Coord _array|count] + point count); 
contro] matrix{count|[0] = Where -> x; 


control matrix(count]{1}] = Where -> y; 
control matrix{count](2] = Where -> z; 


} 
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/** Generate the matrix to compute the forward difference on. 
The forward difference matrix is equal to 
Precision Matrix V xcurrent V_ basis x control matrix. **/ 


pushmatrix (); 

loadmatrix (control matrix); 
multmatrix (current V_ basis): 
multmatrix (Precision Matrix V): 
getmatrix (inter1); 

popmatrix (); 


/* Generate the points on the curve in the V direction */ 
for (t count = 0;t count <= VCURVES:t count++. total points++) { 


(Patch array + total points) -> x = inter1(3]/0]: 
(Patch array + total points) -> y = inter1(3](1]: 
(Patch array + total points) -> z = inter1/[3][2]; 


/* Do the forward difference. * / 
for (Row = 3; Row > 0: Row--) 
for (Column = 0: Column < 4; Column+-+) 
inter1{Row]/Column] = inter1{Row]|[Column] + inter1/Row - 1](Column]; 
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/* Place to put a triangle to send out to user */ 
Coord Triangle 1[4][3]; 


register Point * Where; 
register unsigned int Row: 
register unsigned int Column: 


/* Decompose the patch into its individual triangles. */ 
for (Row = 0; Row < UCURVES: Row+-+) 
for (Column = 0; Column < VCURVES; Column-+-+) { 


Where = (Patch array + ((VCURVES + 1) * Row) + Column); 


Triangle 1(0][(0] = Where -> x; 
Triangle 1(0][1] = Where -> y: 
Triangle 1(0](2] = Where -> z; 


Triangle 1[1] 
Triangle 1{1] 
Triangle 1{1] 


(Where + 1) - 
(Where + 1) -> y: 
(Where + 1) -> z; 


[0] 
[1] 
[2] 


Triangle 1[2][0] 
Triangle 1{2][1] 
Triangle 1/2]/[2] 


(Where + VCURVES + 1) -> x; 
(Where + VCURVES + 1) -> y; 
(Where + VCURVES + 1) -> z; 


Triangle 1{3][(0] = (Where + VCURVES + 2) -> x; 
Triangle 1{3][1] = (Where + VCURVES + 2) -> y; 
Triangle 1{3][2] = (Where + VCURVES + 2) - 


/* Does the user have an intercept routine? * / 
if (_ User routine is on) { 
/* Yes */ 
(* User routine) (&Triangle 1[0]([0]); 
(* User routine) (&Triangle _1[1][0]); 
} 
else { 
/* No user routine, so we use the default. */ 
Default system routine (&Triangle 1[0]/0]); 
Default system routine (&Triangle _1[1]}[0]); 
} 
} 
} 


/* Return the memory used to the system. */ 
cfree (Patch array, (UCURVES + 1) * (VCURVES + 1), sizeof (Point)); 


} 
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npatch 


This function rearranges the input matrices into a form readily 
used in computing points along the four curves defined by those 
matrices. It then computes points for each curve in the U 
direction using the technique of forwards differencing of a 
matrix. Using these points we can then generate points along 
the a curve in the V direction. 
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void npatch (geomx, geomy, geomz) 
Coord geomx/(4]/[4], geomy[4][4], geomz[4][4]; 


register Point * Coord array|4]; 


/* One control matrix for each curve. */ 
Matrix ctr] pts]; 
Matrix ctrl pts2; 
Matrix ctrl pts3; 
Matrix ctrl pts4; 
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/* Load the appropriate control matrix for each curve. */ 


/* Curve 1 */ 

ctrl pts1(0][0] = geomx[0][0]; 
ctr] pts1[0][1] = geomy(0][0]; 
ctrl pts1[0][2] = geomz[0][0]; 


ctrl pts1[{1][0] = geomx[1][0]; 
ctrl pts1[1][1] = geomy(1][0); 
ctrl pts1[{1][2] = geomz(1][0): 


ctrl pts1[2][0] = geomx[2][0); 
ctrl pts1[2][1] = geomy([2](0]; 
ctrl pts1[2][2] = geomz[2][0}: 


ctrl pts1[3](0] = geomx[3][0}; 
ctrl pts1[3][1] = geomy([3] (0): 
ctrl pts1[3][2] = geomz(3][0]: 


/* Curve 2 */ 

ctrl pts2[0](0] = geomx/(0](1); 
ctrl pts2[0][1] = geomy(0][1); 
ctrl pts2(0](2] = geomz([0][1); 


ctrl pts2[1][0] = geomx{1][1]; 
ctrl pts2[1][1] = geomy(1](1]; 
ctrl pts2[1][2] = geomz[1][1]; 


ctrl pts2[2][0] = geomx(2][1); 
ctrl pts2/2](1] = geomy(2][1); 
ctrl pts2[2][2] = geomz(2][1)}; 


ctrl pts2[3][(0] = geomx([3][1]; 


ctrl pts2[3][1] = geomy/[3][1]; 
ctrl pts2[3][2] = geomz(3][1); 
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/* Curve 3 */ 


ctrl pts3[0][0] = geomx(0][2}: 
ctrl pts3[0][1] = geomy|0][2]: 


ctr] pts3[0][2] = geomz[0}[2]: 


ctr] pts3[1][0] = geomx[1][2]; 
ctr] pts3[1][1] = geomy(1][2]: 


ctrl] pts3[1][2] = geomz[1][2}]: 


ctrl pts3[2][0] = geomx[2][2]; 
ctr] pts3[2][1] = geomy[2]|[2]; 


ctrl pts3[2][2] = geomz[2][2]: 


ctr] pts3[3][0] = geomx[3][2]: 
ctr] pts3[3][1] = geomy([3][2]; 


ctrl pts3[3][2] = geomz[3][2]: 


/* Curve 4 */ 


ctr] pts4[0][0] = geomx|[0][3]; 
ctrl pts4/0][1] = geomy/(0][3]; 


ctrl pts4[0][2] = geomz[0][3}: 


ctr] pts4[1][0] = geomx|1][3]; 
ctrl pts4[1][1] = geomy|1]([3]; 


ctrl pts4[1][2] = geomz[1][3]; 


ctr] pts4[2][0] = geomx|2][3]; 
ctrl pts4[2][1] = geomy([2][3]; 


ctr] pts4[2][2] = geomz/[2][3]: 


ctr] pts4[3][0] = geomx[3][3}; 
ctr] pts4[3][1] = geomy([3][3]; 


ctr] pts4[3][2] = geomz[3][3]; 
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register Point *Where; 


register unsigned int Row; 
register unsigned int Column; 
register unsigned int point count; 
register unsigned int count; 


/* An array of pointers to our control matrices. * / 
fioat *matrix pointer|4]; 


/* Matrix used to hold mathematical results. * / 
Matrix inter1l; 


/* Initialize the array */ 
matrix pointer(0] = (float *) ctrl pts1; 


matrix pointer|1] = (float *) ctrl pts2; 
matrix pointer|2] = (float *) ctrl pts3; 
matrix pointer[3] = (float *) ctrl pts4: 


/* Get enough memory to hold the points. */ 
Coord array|[0] = (Point *) calloc (UCURVES + 1. sizeof (Point)): 
Coord array[1] = (Point *) calloc (UCURVES + 1, sizeof (Point)):; 
Coord array[2] = (Point *) calloc (UCURVES + 1, sizeof (Point)): 
Coord array[3] = (Point *) calloc (UCURVES + 1, sizeof (Point)); 
/* For each curve. */ 
for (count = 0; count < 4; count+-+) { 


/* Generate the matrix used in the forward difference for this curve. */ 
pushmatrix (); 
loadmatrix (matrix pointer[count]); 
multmatrix (current U_ basis); 
multmatrix (Precision Matrix U); 
getmatrix (inter1); 
popmatrix (); 
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/* For each curve generate points in the U parametric direction. */ 
for (point count = 0; point count <= UCURVES; point count+-+) { 


Where = (Point *) (Coord array|count] + point count): 


Where -> x = inter1[3]/0]; 
Where -> y = inter1(3](1]; 
Where -> z = inter1[3][2]; 


/* Do the forward difference. * / 

for (Row = 3: Row > 0; Row--) 
for (Column = 0; Column < 4; Column+-+) 
inter1{[Row]{Column] = inter1/Row][Column] + interi{Row - 1]/[Column]: 


} 


/* Call function to finish computations and display. */ 
nspeckle (Coord array); 


/* Return the memory used back to the system. */ 
cfree (Coord array|0], (UCURVES + 1), sizeof (Point) 


" ); 
cfree (Coord array(1],(UCURVES + 1), sizeof (Point)); 
cfree (Coord array|2],(UCURVES + 1), sizeof (Point)); 
F 


cfree (Coord array|3],(UCURVES + 1), sizeof (Point) ); 
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